stanEx3
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}ex9
explanatory variable obs.x also has noise
x~N(x0.sx)
y~N(b0+b1*x0,s) ex8-0.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -47080
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 21 -39.0024 0.0023489 0.000547583 1 1 77
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -39.00
## b0 8.41
## b1 1.68
## s 4.26
## m0[1] 9.31
## m0[2] 33.85
## m0[3] 24.79
## m0[4] 38.10
## m0[5] 7.71
## m0[6] 25.79
## m0[7] 5.55
## m0[8] 27.47
## m0[9] 41.37
## m0[10] 36.98
## m0[11] 4.55
## m0[12] 11.87
## m0[13] 33.88
## m0[14] 19.90
## m0[15] 22.36
## m0[16] 34.77
## m0[17] 18.98
## m0[18] 39.96
## m0[19] 29.58
## m0[20] 19.55
## y0[1] 12.18
## y0[2] 37.30
## y0[3] 21.40
## y0[4] 39.16
## y0[5] 9.97
## y0[6] 28.36
## y0[7] 14.26
## y0[8] 27.45
## y0[9] 43.77
## y0[10] 40.46
## y0[11] -0.78
## y0[12] 10.67
## y0[13] 32.46
## y0[14] 23.24
## y0[15] 22.81
## y0[16] 33.88
## y0[17] 20.17
## y0[18] 36.20
## y0[19] 32.55
## y0[20] 19.58
## m1[1] 4.55
## m1[2] 8.65
## m1[3] 12.74
## m1[4] 16.83
## m1[5] 20.92
## m1[6] 25.01
## m1[7] 29.10
## m1[8] 33.19
## m1[9] 37.28
## m1[10] 41.37
## y1[1] -8.70
## y1[2] 8.55
## y1[3] 9.77
## y1[4] 22.88
## y1[5] 20.64
## y1[6] 29.07
## y1[7] 24.51
## y1[8] 30.62
## y1[9] 36.82
## y1[10] 44.40
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -39.19 -38.84 1.37 1.13 -41.88 -37.72 1.02 388 541
## b0 8.51 8.52 1.99 1.91 5.28 11.72 1.02 201 210
## b1 1.67 1.67 0.17 0.16 1.40 1.95 1.02 304 420
## s 4.83 4.69 0.89 0.79 3.65 6.43 1.01 1100 1088
## m0[1] 9.40 9.42 1.92 1.84 6.32 12.50 1.02 201 212
## m0[2] 33.81 33.84 1.41 1.34 31.44 36.07 1.01 2420 1450
## m0[3] 24.80 24.83 1.10 1.09 22.99 26.53 1.01 513 698
## m0[4] 38.04 38.07 1.70 1.61 35.24 40.78 1.01 1618 1362
## m0[5] 7.81 7.83 2.05 1.97 4.45 11.12 1.02 202 212
## m0[6] 25.79 25.80 1.10 1.09 23.97 27.51 1.01 658 864
## m0[7] 5.66 5.68 2.23 2.17 1.96 9.30 1.02 204 220
## m0[8] 27.47 27.49 1.12 1.12 25.59 29.21 1.01 1008 1276
## m0[9] 41.29 41.33 1.97 1.86 38.00 44.45 1.01 1131 1325
## m0[10] 36.93 36.95 1.62 1.53 34.27 39.54 1.01 1806 1327
## m0[11] 4.68 4.69 2.32 2.26 0.80 8.45 1.02 204 222
## m0[12] 11.95 11.96 1.71 1.67 9.17 14.74 1.02 203 233
## m0[13] 33.84 33.86 1.41 1.34 31.46 36.10 1.01 2414 1450
## m0[14] 19.94 19.95 1.21 1.14 17.90 21.91 1.01 261 283
## m0[15] 22.39 22.41 1.13 1.10 20.51 24.20 1.01 329 366
## m0[16] 34.72 34.75 1.47 1.37 32.27 37.08 1.01 2254 1466
## m0[17] 19.03 19.02 1.25 1.20 16.87 21.07 1.01 247 254
## m0[18] 39.89 39.92 1.85 1.75 36.79 42.85 1.01 1304 1344
## m0[19] 29.56 29.58 1.19 1.18 27.57 31.41 1.01 1852 1528
## m0[20] 19.59 19.60 1.23 1.17 17.50 21.59 1.01 255 257
## y0[1] 9.43 9.42 5.13 5.01 0.90 17.89 1.00 917 1560
## y0[2] 33.79 33.86 5.11 4.83 25.15 41.94 1.00 1778 1778
## y0[3] 24.79 24.89 4.93 4.98 16.85 32.80 1.00 1385 1464
## y0[4] 38.14 38.16 5.33 4.96 29.21 46.74 1.00 1868 1814
## y0[5] 7.74 7.64 5.18 5.16 -0.44 16.19 1.00 854 1592
## y0[6] 25.83 25.94 5.11 5.06 17.70 34.01 1.00 2015 1888
## y0[7] 5.76 5.68 5.49 5.24 -3.47 14.79 1.00 893 1354
## y0[8] 27.28 27.26 5.22 4.83 18.83 35.69 1.00 1786 1965
## y0[9] 41.28 41.26 5.29 5.15 32.55 49.91 1.00 1794 1968
## y0[10] 36.91 36.74 5.15 5.05 28.63 45.51 1.00 1941 1845
## y0[11] 4.74 4.61 5.45 5.17 -3.93 13.83 1.01 674 1579
## y0[12] 12.07 12.12 5.28 5.09 3.29 20.76 1.00 1192 1727
## y0[13] 33.81 33.72 5.31 5.15 25.16 42.30 1.00 1906 1818
## y0[14] 19.89 19.82 5.00 4.67 11.77 28.39 1.00 1783 1973
## y0[15] 22.24 22.14 5.03 4.86 14.01 30.36 1.00 2022 1892
## y0[16] 34.70 34.76 5.10 4.73 26.33 42.68 1.00 2133 1938
## y0[17] 18.94 19.09 4.89 4.67 10.73 26.71 1.00 1503 1482
## y0[18] 39.96 39.94 5.13 4.93 31.48 48.06 1.00 1918 2058
## y0[19] 29.57 29.61 5.00 4.81 21.52 37.88 1.00 1942 1703
## y0[20] 19.62 19.55 5.11 4.95 11.33 27.94 1.00 1947 1780
## m1[1] 4.68 4.69 2.32 2.26 0.80 8.45 1.02 204 222
## m1[2] 8.74 8.76 1.97 1.88 5.55 11.92 1.02 201 210
## m1[3] 12.81 12.81 1.65 1.61 10.13 15.50 1.02 205 236
## m1[4] 16.88 16.87 1.37 1.33 14.58 19.12 1.02 224 244
## m1[5] 20.95 20.97 1.17 1.12 19.00 22.85 1.01 283 306
## m1[6] 25.02 25.03 1.10 1.08 23.20 26.74 1.01 544 706
## m1[7] 29.08 29.10 1.17 1.17 27.09 30.90 1.01 1590 1524
## m1[8] 33.15 33.18 1.37 1.31 30.83 35.33 1.01 2521 1513
## m1[9] 37.22 37.24 1.64 1.54 34.53 39.86 1.01 1753 1346
## m1[10] 41.29 41.33 1.97 1.86 38.00 44.45 1.01 1131 1325
## y1[1] 4.78 4.60 5.54 5.19 -4.41 14.23 1.00 942 1084
## y1[2] 8.84 8.86 5.40 5.16 0.13 17.50 1.00 874 1663
## y1[3] 12.71 12.75 5.23 5.15 4.42 21.28 1.00 849 1331
## y1[4] 16.71 16.73 5.15 4.91 8.21 24.97 1.00 1559 1695
## y1[5] 21.01 20.96 5.08 4.81 12.87 29.54 1.00 1883 1950
## y1[6] 25.03 24.99 4.95 4.69 16.87 33.21 1.00 1852 1712
## y1[7] 29.11 29.12 5.10 4.67 20.58 37.21 1.00 2023 1872
## y1[8] 33.06 33.02 5.20 4.99 24.44 41.33 1.00 2048 1945
## y1[9] 37.05 37.05 5.11 5.05 28.91 45.23 1.00 1907 1786
## y1[10] 41.43 41.27 5.42 5.21 32.61 50.18 1.00 2118 1801
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -45.00 -48.12 13.28 11.60 -62.44 -19.59 1.04 82 61
## b0 8.47 8.40 1.92 1.72 5.29 11.66 1.01 566 968
## b1 1.67 1.68 0.17 0.14 1.41 1.96 1.01 407 937
## s 3.32 3.50 1.69 1.84 0.49 6.01 1.08 41 38
## sx 1.86 1.90 1.11 1.23 0.27 3.54 1.04 73 107
## x0[1] 2.11 1.94 1.70 1.99 -0.21 4.78 1.03 105 1193
## x0[2] 18.52 18.18 2.97 3.87 14.85 23.10 1.04 86 375
## x0[3] 9.00 9.09 1.27 1.18 6.89 10.88 1.01 353 1752
## x0[4] 18.57 18.39 1.52 1.24 16.53 21.21 1.01 286 545
## x0[5] -1.69 -1.38 1.76 1.58 -4.75 0.52 1.02 222 609
## x0[6] 10.40 10.41 1.18 0.83 8.41 12.33 1.01 2270 1360
## x0[7] -0.97 -1.17 1.52 1.28 -3.08 1.40 1.02 227 953
## x0[8] 10.33 10.45 1.39 1.41 8.03 12.39 1.02 172 977
## x0[9] 19.37 19.38 1.36 1.03 17.24 21.60 1.01 521 682
## x0[10] 16.66 16.71 1.29 0.96 14.67 18.70 1.01 701 1375
## x0[11] -2.08 -2.10 1.42 1.10 -4.38 0.11 1.01 410 1037
## x0[12] 3.46 3.36 1.63 1.94 1.17 6.00 1.03 132 1493
## x0[13] 14.70 14.72 1.28 1.07 12.66 16.73 1.01 428 1525
## x0[14] 5.36 5.47 1.66 1.85 2.62 7.60 1.03 123 1060
## x0[15] 8.62 8.59 1.13 0.94 6.77 10.50 1.00 767 1193
## x0[16] 16.14 15.96 1.32 0.98 14.24 18.48 1.01 601 491
## x0[17] 3.89 4.10 2.20 2.82 0.29 6.71 1.04 100 717
## x0[18] 18.76 18.71 1.35 0.94 16.69 21.03 1.01 1083 856
## x0[19] 11.94 12.03 1.35 1.19 9.81 13.91 1.02 245 2039
## x0[20] 6.39 6.43 1.15 0.83 4.44 8.22 1.02 1467 1650
## m0[1] 12.05 12.25 2.90 3.49 7.27 16.15 1.03 121 984
## m0[2] 39.32 39.06 4.65 6.40 32.61 45.86 1.05 74 712
## m0[3] 23.53 23.47 2.09 2.07 20.31 26.91 1.01 316 2501
## m0[4] 39.43 39.64 2.30 2.24 35.51 42.81 1.01 343 1837
## m0[5] 5.75 5.46 2.70 2.85 1.96 10.26 1.02 188 1365
## m0[6] 25.86 25.83 1.93 1.51 22.84 28.99 1.02 3717 1838
## m0[7] 6.92 7.15 2.50 2.31 2.51 10.75 1.01 325 1142
## m0[8] 25.75 25.69 2.31 2.48 22.24 29.52 1.02 184 2505
## m0[9] 40.77 40.65 2.25 1.86 37.06 44.47 1.02 1536 1814
## m0[10] 36.28 36.14 2.14 1.87 33.00 39.82 1.01 1141 2220
## m0[11] 5.10 5.27 2.32 1.88 1.09 8.82 1.02 1055 1477
## m0[12] 14.29 14.54 2.78 3.11 9.67 18.24 1.03 137 1275
## m0[13] 33.00 32.87 2.09 1.91 29.77 36.53 1.01 635 2381
## m0[14] 17.49 17.52 2.64 3.22 13.70 21.67 1.03 119 2546
## m0[15] 22.89 22.97 1.90 1.66 19.82 25.91 1.01 1286 2022
## m0[16] 35.40 35.48 2.08 1.82 31.86 38.69 1.00 1025 2344
## m0[17] 15.05 15.11 3.54 4.71 10.05 20.39 1.04 90 1443
## m0[18] 39.76 39.84 2.18 1.77 36.11 43.33 1.01 2392 2274
## m0[19] 28.41 28.30 2.17 2.04 25.15 32.00 1.02 266 2527
## m0[20] 19.19 19.08 1.87 1.56 16.26 22.34 1.01 2294 2455
## y0[1] 12.11 12.95 4.64 4.10 3.66 18.34 1.01 416 2034
## y0[2] 39.39 40.19 5.82 6.57 28.84 46.79 1.03 124 1128
## y0[3] 23.50 23.03 4.19 3.30 17.17 31.03 1.01 2515 3194
## y0[4] 39.43 40.00 4.38 3.50 31.74 46.15 1.02 1601 1294
## y0[5] 5.74 4.94 4.62 3.78 -0.94 13.80 1.01 1190 2883
## y0[6] 25.91 25.91 4.24 3.18 18.86 33.27 1.02 4117 2972
## y0[7] 7.02 7.52 4.45 3.56 -0.78 13.98 1.02 1609 2660
## y0[8] 25.74 25.11 4.30 3.49 19.45 33.44 1.01 920 2042
## y0[9] 40.88 40.62 4.33 3.42 33.92 48.28 1.02 3250 2678
## y0[10] 36.32 36.03 4.18 3.36 29.60 43.56 1.02 3885 2613
## y0[11] 5.13 5.37 4.35 3.29 -2.35 12.34 1.03 2857 2069
## y0[12] 14.29 15.17 4.62 3.74 5.67 20.57 1.01 404 1851
## y0[13] 32.95 32.62 4.25 3.37 26.31 40.33 1.02 2656 3015
## y0[14] 17.48 16.61 4.57 3.74 11.18 25.96 1.01 672 1851
## y0[15] 22.84 23.19 4.24 3.29 15.36 29.63 1.02 3672 3128
## y0[16] 35.40 35.70 4.27 3.37 27.91 42.15 1.02 3333 3434
## y0[17] 15.01 14.10 5.12 5.03 8.32 24.31 1.02 191 1746
## y0[18] 39.67 39.73 4.33 3.42 32.37 46.61 1.03 3156 2625
## y0[19] 28.46 27.89 4.35 3.41 21.87 36.18 1.02 1576 2671
## y0[20] 19.23 18.99 4.17 3.18 12.54 26.39 1.03 3679 2730
## m1[1] 4.65 4.56 2.25 1.98 0.91 8.34 1.00 527 959
## m1[2] 8.71 8.63 1.90 1.71 5.58 11.86 1.01 569 968
## m1[3] 12.78 12.72 1.58 1.38 10.20 15.38 1.01 653 901
## m1[4] 16.84 16.77 1.31 1.15 14.71 19.04 1.01 833 1150
## m1[5] 20.90 20.82 1.13 0.99 19.10 22.79 1.01 1098 1314
## m1[6] 24.97 24.92 1.06 0.95 23.26 26.74 1.01 1364 1424
## m1[7] 29.03 29.01 1.15 1.05 27.21 30.96 1.01 1205 1454
## m1[8] 33.10 33.13 1.35 1.21 30.90 35.26 1.02 937 1087
## m1[9] 37.16 37.26 1.63 1.45 34.45 39.76 1.02 662 902
## m1[10] 41.22 41.36 1.96 1.70 37.95 44.33 1.02 542 880
## x1[1] -2.27 -2.26 2.16 1.40 -5.75 1.26 1.02 4003 2466
## x1[2] 0.13 0.15 2.12 1.42 -3.36 3.72 1.02 3776 1972
## x1[3] 2.62 2.58 2.14 1.45 -0.84 6.24 1.02 4080 1881
## x1[4] 4.96 4.98 2.17 1.36 1.35 8.63 1.02 3985 2895
## x1[5] 7.42 7.42 2.21 1.44 3.82 11.01 1.02 3914 1734
## x1[6] 9.84 9.85 2.15 1.39 6.37 13.41 1.02 3945 1714
## x1[7] 12.27 12.29 2.15 1.44 8.68 15.85 1.02 4218 1823
## x1[8] 14.76 14.72 2.15 1.44 11.36 18.41 1.01 3977 2317
## x1[9] 17.14 17.18 2.23 1.41 13.57 20.77 1.02 4113 1912
## x1[10] 19.60 19.59 2.22 1.47 15.83 23.22 1.02 3867 2322
## y1[1] 4.63 4.42 4.40 3.64 -2.24 12.08 1.02 1642 2065
## y1[2] 8.68 8.46 4.08 3.40 2.08 15.48 1.02 1730 2227
## y1[3] 12.81 12.71 4.11 3.33 6.13 19.53 1.01 1907 2417
## y1[4] 16.85 16.74 3.94 3.14 10.42 23.25 1.02 3393 2750
## y1[5] 20.94 20.86 3.88 3.06 14.30 27.23 1.02 3606 2506
## y1[6] 24.91 24.88 3.78 2.96 18.47 31.23 1.02 3973 3161
## y1[7] 29.00 29.07 3.87 3.04 22.40 35.30 1.03 3454 1617
## y1[8] 33.09 33.05 4.04 3.16 26.48 39.72 1.02 2508 2743
## y1[9] 37.21 37.29 4.04 3.29 30.53 43.93 1.02 2510 2734
## y1[10] 41.11 41.30 4.22 3.51 34.04 47.98 1.02 1593 2332
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
ex10
outlier
y~cauchy(b0+b1*x,s)ex10.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -1673.04
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -32.8956 0.0011713 0.000402703 0.8809 0.8809 29
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -32.90
## b0 7.11
## b1 1.67
## s 3.14
## m0[1] 12.13
## m0[2] 20.03
## m0[3] 12.76
## m0[4] 9.63
## m0[5] 10.58
## m0[6] 16.27
## m0[7] 10.69
## m0[8] 19.16
## m0[9] 17.03
## m0[10] 8.85
## m0[11] 18.34
## m0[12] 19.89
## m0[13] 13.73
## m0[14] 16.61
## m0[15] 13.81
## m0[16] 16.51
## m0[17] 9.57
## m0[18] 15.16
## m0[19] 22.00
## m0[20] 21.08
## y0[1] 10.47
## y0[2] 18.10
## y0[3] 14.27
## y0[4] 12.60
## y0[5] 11.86
## y0[6] 16.84
## y0[7] 8.54
## y0[8] 20.21
## y0[9] 13.24
## y0[10] 8.07
## y0[11] 20.48
## y0[12] 18.28
## y0[13] 10.81
## y0[14] 16.93
## y0[15] 16.76
## y0[16] 13.20
## y0[17] 9.40
## y0[18] 16.45
## y0[19] 24.08
## y0[20] 18.91
## m1[1] 8.85
## m1[2] 10.31
## m1[3] 11.77
## m1[4] 13.23
## m1[5] 14.69
## m1[6] 16.15
## m1[7] 17.61
## m1[8] 19.08
## m1[9] 20.54
## m1[10] 22.00
## y1[1] 6.84
## y1[2] 11.54
## y1[3] 7.90
## y1[4] 11.21
## y1[5] 15.30
## y1[6] 14.52
## y1[7] 17.47
## y1[8] 19.48
## y1[9] 18.87
## y1[10] 23.52
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.29 -32.98 1.22 1.02 -35.54 -31.95 1.01 662 878
## b0 7.23 7.22 1.66 1.72 4.50 9.97 1.02 448 411
## b1 1.65 1.66 0.32 0.32 1.15 2.18 1.01 515 543
## s 3.53 3.46 0.64 0.60 2.65 4.69 1.00 1099 1340
## m0[1] 12.18 12.19 0.94 0.96 10.58 13.72 1.01 595 796
## m0[2] 19.99 19.99 1.25 1.17 17.92 22.04 1.00 1317 1450
## m0[3] 12.81 12.81 0.88 0.90 11.33 14.25 1.01 677 829
## m0[4] 9.72 9.72 1.26 1.29 7.57 11.81 1.01 472 435
## m0[5] 10.65 10.66 1.13 1.14 8.75 12.52 1.01 498 536
## m0[6] 16.27 16.26 0.83 0.80 14.94 17.65 1.00 2333 1397
## m0[7] 10.77 10.77 1.11 1.13 8.89 12.60 1.01 503 555
## m0[8] 19.12 19.12 1.13 1.07 17.27 21.00 1.00 1570 1428
## m0[9] 17.02 17.02 0.89 0.86 15.58 18.47 1.00 2427 1308
## m0[10] 8.94 8.94 1.38 1.41 6.64 11.25 1.02 459 429
## m0[11] 18.32 18.31 1.02 0.98 16.64 19.98 1.00 1991 1421
## m0[12] 19.85 19.85 1.23 1.16 17.83 21.86 1.00 1350 1464
## m0[13] 13.76 13.76 0.82 0.83 12.39 15.10 1.00 884 1138
## m0[14] 16.61 16.60 0.85 0.81 15.24 18.02 1.00 2401 1401
## m0[15] 13.84 13.84 0.81 0.82 12.47 15.17 1.00 907 1151
## m0[16] 16.51 16.51 0.85 0.80 15.15 17.91 1.00 2385 1422
## m0[17] 9.65 9.65 1.27 1.30 7.49 11.76 1.01 471 433
## m0[18] 15.17 15.18 0.79 0.77 13.85 16.47 1.00 1592 1384
## m0[19] 21.92 21.92 1.56 1.50 19.31 24.48 1.01 1009 1407
## m0[20] 21.02 21.01 1.41 1.34 18.65 23.29 1.00 1122 1471
## y0[1] 12.32 12.25 3.71 3.55 6.33 18.43 1.00 1994 2000
## y0[2] 20.04 20.05 3.86 3.63 13.69 26.40 1.00 1817 1805
## y0[3] 12.90 12.88 3.64 3.51 6.95 18.60 1.00 1672 1930
## y0[4] 9.69 9.76 3.82 3.69 3.35 15.73 1.00 1556 1829
## y0[5] 10.56 10.60 3.79 3.69 4.16 16.68 1.00 1518 1705
## y0[6] 16.17 16.20 3.72 3.47 10.21 22.22 1.00 2090 1877
## y0[7] 10.66 10.67 3.67 3.63 4.87 16.67 1.00 1552 1872
## y0[8] 19.23 19.23 3.69 3.59 13.21 25.24 1.00 2002 2042
## y0[9] 17.00 17.09 3.64 3.45 11.08 22.92 1.00 2116 1926
## y0[10] 8.82 8.91 3.79 3.50 2.66 15.04 1.00 1390 1916
## y0[11] 18.20 18.24 3.76 3.70 12.05 24.34 1.00 2107 1969
## y0[12] 19.63 19.74 3.82 3.74 13.26 25.84 1.00 2066 1848
## y0[13] 13.81 13.87 3.70 3.74 7.74 19.93 1.00 1947 1864
## y0[14] 16.65 16.76 3.74 3.74 10.39 22.66 1.00 1892 1794
## y0[15] 13.71 13.74 3.81 3.62 7.50 19.94 1.00 1891 2080
## y0[16] 16.69 16.66 3.68 3.58 10.55 22.83 1.00 1967 1964
## y0[17] 9.56 9.57 3.88 3.64 3.14 15.95 1.00 1283 1313
## y0[18] 15.22 15.29 3.71 3.56 8.87 21.38 1.00 2147 1954
## y0[19] 21.85 21.93 3.99 3.65 15.27 28.31 1.00 1728 1651
## y0[20] 21.04 21.08 3.91 3.66 14.55 27.33 1.00 1524 1848
## m1[1] 8.94 8.94 1.38 1.41 6.64 11.25 1.02 459 429
## m1[2] 10.39 10.39 1.17 1.18 8.41 12.32 1.01 489 507
## m1[3] 11.83 11.85 0.98 0.99 10.17 13.42 1.01 563 692
## m1[4] 13.27 13.27 0.85 0.86 11.85 14.65 1.01 760 1053
## m1[5] 14.71 14.72 0.79 0.79 13.37 16.00 1.00 1296 1307
## m1[6] 16.16 16.15 0.82 0.80 14.82 17.51 1.00 2297 1439
## m1[7] 17.60 17.59 0.94 0.91 16.04 19.15 1.00 2381 1243
## m1[8] 19.04 19.05 1.12 1.05 17.20 20.90 1.00 1604 1428
## m1[9] 20.48 20.49 1.33 1.25 18.31 22.64 1.00 1214 1449
## m1[10] 21.92 21.92 1.56 1.50 19.31 24.48 1.01 1009 1407
## y1[1] 8.93 8.92 3.83 3.80 2.65 15.39 1.00 1547 1729
## y1[2] 10.31 10.34 3.72 3.69 4.27 16.32 1.00 1613 1785
## y1[3] 11.83 11.86 3.59 3.37 5.90 17.78 1.00 1616 1738
## y1[4] 13.25 13.29 3.71 3.60 7.07 19.51 1.00 1884 1954
## y1[5] 14.65 14.71 3.76 3.69 8.34 20.73 1.00 1719 1914
## y1[6] 16.13 16.02 3.69 3.61 10.08 21.93 1.00 1899 2013
## y1[7] 17.62 17.55 3.72 3.39 11.73 23.78 1.00 1987 2054
## y1[8] 19.19 19.22 3.79 3.73 12.76 25.17 1.00 1944 1934
## y1[9] 20.37 20.38 3.93 3.76 13.96 26.95 1.00 1842 1853
## y1[10] 21.94 21.94 3.96 3.90 15.52 28.56 1.00 1918 1620
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -118.839
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -15.0002 7.45559e-05 0.000492157 0.8905 0.8905 35
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -15.00
## b0 5.69
## b1 1.84
## s 0.75
## m0[1] 11.20
## m0[2] 19.89
## m0[3] 11.90
## m0[4] 8.46
## m0[5] 9.50
## m0[6] 15.75
## m0[7] 9.62
## m0[8] 18.92
## m0[9] 16.59
## m0[10] 7.59
## m0[11] 18.03
## m0[12] 19.73
## m0[13] 12.96
## m0[14] 16.12
## m0[15] 13.05
## m0[16] 16.01
## m0[17] 8.39
## m0[18] 14.53
## m0[19] 22.04
## m0[20] 21.03
## y0[1] 11.85
## y0[2] 20.25
## y0[3] 10.37
## y0[4] 7.82
## y0[5] 9.65
## y0[6] 16.38
## y0[7] 9.46
## y0[8] 18.58
## y0[9] 12.31
## y0[10] 10.16
## y0[11] -44.34
## y0[12] 19.90
## y0[13] 13.45
## y0[14] 15.79
## y0[15] 15.72
## y0[16] 16.37
## y0[17] 7.54
## y0[18] 15.08
## y0[19] 16.34
## y0[20] 20.49
## m1[1] 7.59
## m1[2] 9.20
## m1[3] 10.81
## m1[4] 12.41
## m1[5] 14.02
## m1[6] 15.62
## m1[7] 17.23
## m1[8] 18.83
## m1[9] 20.44
## m1[10] 22.04
## y1[1] 8.11
## y1[2] 8.26
## y1[3] 11.13
## y1[4] 13.07
## y1[5] 12.20
## y1[6] 13.64
## y1[7] 15.91
## y1[8] 18.57
## y1[9] 20.32
## y1[10] 23.13
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.89 -16.60 1.29 1.08 -19.41 -15.45 1.00 579 688
## b0 5.52 5.57 0.66 0.65 4.40 6.54 1.02 318 422
## b1 1.85 1.85 0.11 0.10 1.68 2.03 1.01 377 515
## s 0.87 0.84 0.26 0.24 0.51 1.36 1.00 854 677
## m0[1] 11.08 11.11 0.42 0.41 10.35 11.73 1.01 416 693
## m0[2] 19.85 19.85 0.41 0.40 19.17 20.50 1.00 2206 1597
## m0[3] 11.79 11.81 0.39 0.39 11.10 12.41 1.01 461 776
## m0[4] 8.31 8.35 0.53 0.53 7.39 9.12 1.02 335 550
## m0[5] 9.36 9.40 0.48 0.49 8.51 10.11 1.01 355 606
## m0[6] 15.68 15.69 0.33 0.35 15.13 16.20 1.00 1369 1417
## m0[7] 9.49 9.52 0.48 0.48 8.65 10.23 1.01 358 606
## m0[8] 18.87 18.87 0.38 0.38 18.25 19.48 1.00 2257 1508
## m0[9] 16.52 16.52 0.34 0.36 15.96 17.05 1.00 1801 1471
## m0[10] 7.44 7.49 0.57 0.56 6.47 8.32 1.02 327 511
## m0[11] 17.97 17.98 0.36 0.37 17.38 18.53 1.00 2219 1515
## m0[12] 19.69 19.69 0.41 0.40 19.02 20.34 1.00 2220 1546
## m0[13] 12.86 12.88 0.36 0.36 12.22 13.43 1.01 583 842
## m0[14] 16.05 16.06 0.33 0.35 15.51 16.57 1.00 1595 1471
## m0[15] 12.94 12.97 0.36 0.36 12.32 13.51 1.01 597 828
## m0[16] 15.94 15.95 0.33 0.35 15.39 16.46 1.00 1538 1481
## m0[17] 8.24 8.28 0.53 0.53 7.32 9.06 1.02 334 550
## m0[18] 14.44 14.46 0.34 0.35 13.87 14.97 1.01 863 1070
## m0[19] 22.02 22.02 0.50 0.47 21.21 22.82 1.00 1504 1415
## m0[20] 21.01 21.00 0.45 0.44 20.25 21.72 1.00 1934 1411
## y0[1] 12.07 11.13 25.14 1.38 5.12 17.02 1.00 1918 1997
## y0[2] 19.63 19.87 25.01 1.37 14.32 25.78 1.00 1903 1925
## y0[3] 12.81 11.87 43.76 1.41 6.18 17.82 1.00 1790 1930
## y0[4] 8.26 8.36 47.67 1.43 3.28 14.11 1.00 1706 1968
## y0[5] 8.60 9.36 18.96 1.45 3.45 15.00 1.00 1693 1973
## y0[6] 16.74 15.70 23.48 1.38 11.27 21.19 1.00 1946 2055
## y0[7] 9.37 9.42 26.47 1.34 3.43 14.33 1.00 1674 1932
## y0[8] 19.44 18.90 18.09 1.32 13.58 24.09 1.00 1922 1956
## y0[9] 16.93 16.51 55.12 1.42 11.13 21.63 1.00 2021 1843
## y0[10] 6.93 7.45 26.45 1.41 2.11 12.69 1.00 1684 1915
## y0[11] 19.36 17.98 45.64 1.39 12.64 23.34 1.00 2039 1959
## y0[12] 17.93 19.72 37.46 1.36 14.16 24.59 1.00 1888 1931
## y0[13] 13.27 12.86 15.69 1.33 8.12 17.88 1.00 2102 2039
## y0[14] 16.43 15.99 43.79 1.30 10.74 21.12 1.00 1852 1832
## y0[15] 9.27 12.99 214.56 1.37 7.57 19.16 1.00 2095 1971
## y0[16] 16.62 15.94 39.23 1.30 11.36 20.89 1.00 1949 1973
## y0[17] 8.67 8.32 20.13 1.41 4.10 14.63 1.00 1693 2027
## y0[18] 18.90 14.44 163.39 1.32 9.33 20.33 1.00 1905 1850
## y0[19] 20.86 22.06 41.15 1.44 17.02 27.37 1.00 2052 2017
## y0[20] 21.34 21.01 29.23 1.36 16.17 26.20 1.00 2016 1845
## m1[1] 7.44 7.49 0.57 0.56 6.47 8.32 1.02 327 511
## m1[2] 9.06 9.10 0.49 0.50 8.19 9.82 1.01 348 585
## m1[3] 10.68 10.72 0.43 0.43 9.93 11.35 1.01 398 703
## m1[4] 12.30 12.33 0.38 0.38 11.64 12.90 1.01 510 832
## m1[5] 13.92 13.95 0.34 0.35 13.34 14.46 1.01 751 961
## m1[6] 15.54 15.55 0.33 0.35 14.99 16.07 1.00 1295 1352
## m1[7] 17.16 17.17 0.34 0.36 16.61 17.71 1.00 2081 1539
## m1[8] 18.78 18.79 0.38 0.38 18.16 19.39 1.00 2256 1508
## m1[9] 20.40 20.40 0.43 0.41 19.70 21.08 1.00 2139 1455
## m1[10] 22.02 22.02 0.50 0.47 21.21 22.82 1.00 1504 1415
## y1[1] 5.08 7.49 93.65 1.46 1.28 13.51 1.00 1901 1944
## y1[2] 9.41 9.08 30.62 1.45 2.96 13.94 1.00 1737 2056
## y1[3] 11.42 10.74 34.16 1.40 4.86 16.42 1.00 1984 1766
## y1[4] 9.79 12.28 75.51 1.51 6.10 18.30 1.00 1895 2010
## y1[5] 11.54 13.90 160.53 1.38 8.19 19.70 1.00 2033 2040
## y1[6] 16.05 15.58 42.81 1.34 10.17 21.27 1.00 1874 1918
## y1[7] 16.13 17.18 49.36 1.31 11.42 23.08 1.00 1873 1914
## y1[8] 20.54 18.75 67.97 1.40 12.61 24.85 1.00 1854 1953
## y1[9] 20.92 20.44 50.03 1.40 14.60 26.21 1.00 1841 1934
## y1[10] 22.58 22.00 23.62 1.38 15.86 27.32 1.00 1904 1971
ex11
censored
y i=1-N, y~N(m,s)
acutualy ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitiondata=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -5182.13
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -9.04902 0.000449008 1.46199e-05 0.9943 0.9943 44
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.05
## b0 12.41
## b1 2.19
## s 1.66
## m0[1] 19.49
## m0[2] 20.64
## m0[3] 19.32
## m0[4] 16.70
## m0[5] 17.44
## m0[6] 16.77
## m0[7] 18.42
## m0[8] 26.00
## m0[9] 16.65
## y0[1] 20.39
## y0[2] 19.66
## y0[3] 19.25
## y0[4] 16.92
## y0[5] 16.74
## y0[6] 16.10
## y0[7] 18.49
## y0[8] 24.81
## y0[9] 16.58
## m1[1] 12.73
## m1[2] 14.85
## m1[3] 16.96
## m1[4] 19.07
## m1[5] 21.19
## m1[6] 23.30
## m1[7] 25.42
## m1[8] 27.53
## m1[9] 29.64
## m1[10] 31.76
## y1[1] 11.41
## y1[2] 15.89
## y1[3] 14.11
## y1[4] 17.12
## y1[5] 20.59
## y1[6] 22.58
## y1[7] 25.16
## y1[8] 28.41
## y1[9] 31.86
## y1[10] 31.81mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -10.40 -9.98 1.61 1.27 -13.67 -8.72 1.01 364 333
## b0 12.65 12.64 2.16 1.82 9.25 16.04 1.01 311 246
## b1 2.12 2.15 0.65 0.54 1.05 3.14 1.01 338 307
## s 2.32 2.12 0.86 0.64 1.38 3.91 1.01 770 603
## m0[1] 19.51 19.53 0.87 0.72 18.12 20.80 1.00 1318 1050
## m0[2] 20.62 20.63 0.98 0.81 19.04 22.06 1.00 1334 886
## m0[3] 19.34 19.36 0.87 0.72 17.96 20.65 1.00 1229 1025
## m0[4] 16.81 16.80 1.11 0.94 15.07 18.56 1.00 392 334
## m0[5] 17.52 17.52 0.99 0.84 15.99 19.08 1.00 473 476
## m0[6] 16.87 16.87 1.10 0.92 15.16 18.61 1.00 398 349
## m0[7] 18.47 18.48 0.88 0.75 17.07 19.86 1.00 742 762
## m0[8] 25.81 25.92 2.25 1.85 22.17 29.18 1.00 462 518
## m0[9] 16.75 16.75 1.12 0.94 15.01 18.53 1.00 388 333
## y0[1] 19.48 19.55 2.66 2.20 15.32 23.34 1.00 1739 1660
## y0[2] 20.60 20.62 2.71 2.21 16.36 24.76 1.01 1706 1454
## y0[3] 19.30 19.37 2.70 2.22 14.84 23.35 1.00 1941 1428
## y0[4] 16.89 16.90 2.69 2.29 12.61 21.01 1.00 1235 1466
## y0[5] 17.46 17.51 2.68 2.29 13.19 21.79 1.00 1508 1414
## y0[6] 16.81 16.79 2.78 2.38 12.56 21.27 1.00 1089 998
## y0[7] 18.48 18.44 2.63 2.20 14.31 22.92 1.00 1783 1552
## y0[8] 25.81 25.86 3.28 2.83 20.56 30.98 1.00 826 1083
## y0[9] 16.69 16.70 2.66 2.29 12.42 20.93 1.00 1329 1061
## m1[1] 12.96 12.95 2.07 1.76 9.69 16.23 1.01 312 240
## m1[2] 15.01 15.00 1.52 1.29 12.63 17.43 1.00 326 276
## m1[3] 17.06 17.05 1.06 0.89 15.40 18.74 1.00 416 394
## m1[4] 19.10 19.13 0.86 0.72 17.73 20.42 1.00 1089 944
## m1[5] 21.15 21.18 1.07 0.89 19.41 22.76 1.00 1099 817
## m1[6] 23.20 23.28 1.53 1.25 20.68 25.46 1.00 603 580
## m1[7] 25.24 25.34 2.08 1.73 21.86 28.38 1.00 480 518
## m1[8] 27.29 27.40 2.67 2.21 22.94 31.34 1.01 433 465
## m1[9] 29.34 29.49 3.28 2.71 23.94 34.38 1.01 410 470
## m1[10] 31.39 31.56 3.89 3.21 24.99 37.37 1.01 395 452
## y1[1] 13.05 12.94 3.40 2.87 7.77 18.69 1.00 598 704
## y1[2] 15.05 14.96 2.98 2.50 10.58 19.99 1.00 837 1087
## y1[3] 17.13 17.08 2.67 2.41 12.97 21.52 1.00 1151 1168
## y1[4] 19.12 19.10 2.67 2.21 14.98 23.49 1.00 1805 1581
## y1[5] 21.06 21.12 2.75 2.34 16.55 25.36 1.00 1977 1586
## y1[6] 23.23 23.32 2.95 2.57 18.47 27.65 1.00 1208 1152
## y1[7] 25.21 25.20 3.20 2.80 20.20 30.24 1.00 939 1049
## y1[8] 27.30 27.42 3.73 3.21 21.31 32.99 1.00 652 684
## y1[9] 29.34 29.61 3.98 3.45 22.76 35.30 1.00 630 667
## y1[10] 31.30 31.49 4.65 4.13 23.99 38.21 1.00 517 609
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -263.167
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 26 -15.7809 0.000214419 0.000130948 1 1 36
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -15.78
## b0 9.07
## b1 3.14
## s 2.44
## m0[1] 19.22
## m0[2] 20.86
## m0[3] 18.97
## m0[4] 15.23
## m0[5] 16.29
## m0[6] 15.33
## m0[7] 17.69
## m0[8] 28.55
## m0[9] 15.15
## y0[1] 18.42
## y0[2] 17.31
## y0[3] 18.82
## y0[4] 13.22
## y0[5] 15.08
## y0[6] 12.99
## y0[7] 22.07
## y0[8] 30.35
## y0[9] 15.84
## m1[1] 9.54
## m1[2] 12.57
## m1[3] 15.60
## m1[4] 18.63
## m1[5] 21.65
## m1[6] 24.68
## m1[7] 27.71
## m1[8] 30.74
## m1[9] 33.76
## m1[10] 36.79
## y1[1] 14.98
## y1[2] 11.13
## y1[3] 15.60
## y1[4] 12.19
## y1[5] 19.70
## y1[6] 25.30
## y1[7] 30.84
## y1[8] 29.33
## y1[9] 35.88
## y1[10] 30.31mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.59 -16.16 1.40 1.15 -19.41 -15.07 1.00 491 630
## b0 7.86 8.18 2.30 2.06 3.59 11.00 1.01 289 370
## b1 3.52 3.42 0.64 0.56 2.68 4.82 1.01 311 430
## s 3.36 3.08 1.14 0.88 2.04 5.52 1.01 617 572
## m0[1] 19.25 19.27 1.03 0.93 17.54 20.92 1.00 962 861
## m0[2] 21.09 21.04 1.09 0.92 19.38 22.94 1.00 1443 1206
## m0[3] 18.97 18.99 1.03 0.93 17.26 20.59 1.00 881 854
## m0[4] 14.76 14.92 1.31 1.19 12.48 16.62 1.01 365 509
## m0[5] 15.95 16.07 1.19 1.08 13.89 17.66 1.01 424 529
## m0[6] 14.88 15.03 1.29 1.18 12.62 16.71 1.01 369 517
## m0[7] 17.54 17.60 1.07 0.97 15.74 19.13 1.00 586 690
## m0[8] 29.72 29.43 2.18 1.84 26.81 34.01 1.01 582 618
## m0[9] 14.68 14.83 1.32 1.21 12.39 16.53 1.01 362 509
## y0[1] 19.18 19.27 3.69 3.29 13.16 25.05 1.00 1916 1522
## y0[2] 21.16 21.13 3.68 3.25 15.32 27.19 1.00 1968 1716
## y0[3] 18.99 18.97 3.64 3.22 13.19 24.72 1.00 1475 1505
## y0[4] 14.86 15.08 3.78 3.45 8.46 20.55 1.00 1253 1572
## y0[5] 16.01 15.95 3.76 3.21 10.11 21.92 1.00 1455 1412
## y0[6] 14.78 14.90 3.69 3.30 8.81 20.52 1.00 1318 1281
## y0[7] 17.59 17.54 3.83 3.44 11.74 23.42 1.00 1820 1722
## y0[8] 29.82 29.53 4.07 3.60 23.85 37.02 1.00 1136 1101
## y0[9] 14.65 14.79 3.85 3.34 8.38 20.64 1.00 1252 1008
## m1[1] 8.38 8.69 2.21 1.98 4.24 11.40 1.01 291 370
## m1[2] 11.78 11.99 1.69 1.52 8.53 14.14 1.01 307 383
## m1[3] 15.18 15.32 1.26 1.16 12.98 16.98 1.01 382 525
## m1[4] 18.58 18.61 1.04 0.93 16.86 20.18 1.00 786 853
## m1[5] 21.98 21.91 1.15 0.98 20.26 24.02 1.00 1535 1048
## m1[6] 25.38 25.21 1.53 1.30 23.25 28.27 1.00 1070 871
## m1[7] 28.78 28.51 2.03 1.71 26.05 32.69 1.01 636 695
## m1[8] 32.18 31.82 2.58 2.20 28.79 37.24 1.01 492 602
## m1[9] 35.58 35.13 3.16 2.67 31.45 41.93 1.01 430 593
## m1[10] 38.98 38.43 3.75 3.17 34.09 46.48 1.01 397 583
## y1[1] 8.34 8.65 4.38 3.84 0.90 14.74 1.00 798 1102
## y1[2] 11.83 11.93 3.89 3.43 5.19 17.86 1.00 909 1383
## y1[3] 15.08 15.17 3.63 3.32 9.10 20.84 1.00 1421 1560
## y1[4] 18.57 18.69 3.66 3.33 12.35 24.24 1.00 1486 1500
## y1[5] 21.91 21.84 3.71 3.22 16.01 28.04 1.00 1918 1693
## y1[6] 25.50 25.29 4.00 3.43 19.51 32.09 1.00 1634 1590
## y1[7] 28.95 28.61 4.15 3.66 23.02 36.21 1.00 1186 1118
## y1[8] 32.11 31.71 4.41 3.82 25.74 40.00 1.00 1206 1115
## y1[9] 35.64 35.16 4.71 4.20 28.89 43.91 1.00 1018 1254
## y1[10] 38.80 38.17 5.12 4.66 31.58 48.23 1.00 672 825
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
bimodal distribution
mixed normal distributiondata {
int<lower=0> N; // データの数
real y[N]; // 観測データ
}
parameters {
real<lower=0, upper=1> theta; // 分布1と分布2の混合比率
real mu1; // 分布1の平均
real mu2; // 分布2の平均
real<lower=0> sigma1; // 分布1の標準偏差
real<lower=0> sigma2; // 分布2の標準偏差
}
model {
for (n in 1:N) {
target += log_mix(theta,
normal_lpdf(y[n] | mu1, sigma1),
normal_lpdf(y[n] | mu2, sigma2));
}
}y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()
EM algorithm
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -246 100 6 -520 -520
##
## Clustering table:
## 1 2 3
## 29 33 38rst$parameters## $pro
## [1] 0.290 0.331 0.378
##
## $mean
## 1 2 3
## -4.977 0.234 4.842
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 0.925
##
##
## $Vinv
## NULLplot(rst)
